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PACS 71.27.+a - Strongly correlated electron systems; heavy fermions. 

PACS 71 . 10 . Fd - Lattice fermion models 

PACS 71 . 38 . -k - Polarons and electron-phonon interactions 

Abstract. - We introduce a new type of Gutzwiller variational wavefunction for correlated elec- 
trons coupled to phonons, able to treat on equal footing electronic and lattice degrees of freedom. 
We benchmark the wavefunction in the infinite-C/ Hubbard-Hlolstein model away from half-filling 
on a Bethe lattice, where we can directly compare with exact results by Dynamical Mean-Field 
Theory. For this model, we find that variational results agree perfectly well with the exact ones. 
In particular the wavefunction correctly describes the crossover to a heavy polaron gas upon 
increasing the electron-phonon coupling. 
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\ The coherent electron motion is strongly slowed down 
by Coulomb repulsion close to a Mott metal-to-insulator 
transition. In this situation the effects of the electron- 
lattice coupling may be enhanced, since the lattice has 
inore time to adjust following the electron motion. In 
fact, the Mott transition in real systems is often accom- 
panied by substantial structural changes that turn it into 
a strong first-order one. as it happens in the prototypi- 
cal example of V2O3. [1] However, despite its relevance in 
many strongly correlated materials, the electron-phonon 
(e-ph) interaction has been just touched on in the con- 
text of the Mott transition. This has happened partly 
because most theoretical efforts have focused mainly on 
the strong correlation alone, or on the purely e-ph driven 
polaron problem, but also because it is not obvious how 
to deal with the c-ph coupling once the electron-elecron 
repulsion slackens quasiparticle motion so much that the 
phonon dynamics becomes comparable or even faster. For 
this reason, recent attempts to understand the role of the 
e-ph coupling near a Mott transition have mainly limited 
to "exact" numerical simulations using several techniques, 
including Quantum Monte Carlo [2-5], exact diagonaliza- 
tion [6], Density-Matrix Rcnormalization Group [7] and 
Dynamical Mean-Field Theory (DMFT) [8,9] or exten- 
sions [10]. However, besides accurate numerical results. 



it is desirable to have at disposal also analytic or semi- 
analytic approaches, even though approximate, that may 
provide us with simple and intuitive physical insights and 
may be extended to situations where numerical approaches 
become hardly feasible. 

A popular and simple tool for the study of strongly cor- 
related electron systems, so far in the absence of e-ph in- 
teraction, is the Gutzwiller variational approach. In its 
simplest version [11], the Gutzwiller approach amounts 
to variationally modify the relative weights of the local 
electronic configurations of an uncorrelated wavefunction. 
This variational wavefunction is able to describe a Mott 
transition only in the limit of infinite coordination [12], 
where the expectation value of the Hamiltonian can be 
calculated analytically. In finite-coordination lattices, it 
has been recently argued [13] that the Gutzwiller wave- 
function may still describe a Mott transition once a long- 
range Jastrow factor is included. In this paper we propose 
an extension of the Gutzwiller wavefunction to account 
for the c-ph coupling, which we label Gutzwiller Phonon 
Wavefunction (GPW). We test the accuracy of this wave- 
function in the Hubbard-Holstein model on a Bethe lat- 
tice, where analytical results can be obtained and com- 
pared with "exact" DMFT calculations. 
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We consider an Hubbard-Holstcin model 

H = 



(1) 



where c^^ (cj^) and b^^ (^L) a-re annihilation(crcation) op- 
erators at site i for spin-cr electrons and for dispersionless 
phonons of frequency (jJq, respectively, U the on-site Hub- 
bard repulsion while a parameterizes the e-ph coupling. 
The hopping tj ^/z connects nearest neighbor sites on a 
Cayley tree with coordination number z. The calculations 
will actually be performed in the limit z ^ oo, namely 
for the so-called Bethe lattice, which has a semi-circular 
density of states with half-bandwidth D — 2t. The e- 
ph properties are conveniently parameterized through the 
adiabaticity parameter 7 ~ loq/D, and the dimensionless 
coupling A — 2o?ijJqID = 2a^^. The choice of a Bethe 
lattice has the advantage to allow for analytical calcula- 
tions which can be compared directly with DMFT, which 
is exact in infinite-coordination lattices. In the DMFT cal- 
culations we use exact diagonalization as impurity solver, 
with the conduction bath approximated by a finite num- 
ber of levels TVf, = 9. We refer to Refs. [9] for more de- 
tails on the practical implementation of this technique to 
the Hubbard-Holstcin model. The phase diagram of the 
Hamiltonian (1) contains several different phases depend- 
ing on the various parameters A, 7, U jD as well as on 
the electron density. In order to simplify the analysis, we 
will assume a very large, actually infinite, U I D, which ex- 
cludes superconductivity. Moreover we will discard the 
possible occurrence of ferromagnetism very close to half- 
filling. Under these assumptions, we expect, away from 
half-filling, an evolution from a strongly correlated metal 
at small A into a gas of heavy polarons at large A. Previous 
DMFT calculations [9] have shown that this evolution is 
continuous for large (but finite) U , although the crossover 
may be quite sharp, while variational Lang-Firsov trans- 
formations within the slave-boson approach [14] find in- 
stead a first order phase transition between the two metal- 
lic phases [15]. We show that polaron formation occurs as 
a crossover also for U / D — > 00, and that our novel varia- 
tional approach agrees extremely well, also quantitatively 
with exact DMFT, correctly describing the crossover. 

In analogy with purely electronic models [16, 17], we 
introduce a Gutzwillcr wavcfunction that depends also on 
the phonon coordinates Xi ~ {hi + bl)/\/2 through: 
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where [vto) is an uncorrclatcd fcrmionic wavcfunction. 
here assumed to be the non interacting Fermi sea. The 
Gutzwiller correlator is defined by 
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where is the projection operator at site i onto 

the subspace with electron number = 0, 1, 2, P,^ are vari- 
ational parameters while P^*^' are the occupation proba- 
bilities of the Fermi sea, i.e., Pp*"^ = (1 - n-/2)^ Pf^ = 
n(l — n/2), and P'^^ = {n/2)'^, n being the average num- 
ber of electrons per site. Finally, 4>^{xi) are normalized 
functions associated to the different local electronic con- 
figurations. The introduction of this variational freedom is 
the novelty of our approach. We emphasize that the wave- 
function (2) can be easily modified to account for other 
types of c-ph coupling as well as for the phonon disper- 
sion, although that would likely require numerical evalua- 
tion of the energy and optimization of the variational pa- 
rameters. For instance, phonon dispersion can be included 
if 1*0) ^{xi,X2, . . .,xn) 1*0), where ^{xi,X2, • . .,xn) 
is a variational many-body phonon wavcfunction that in- 
cludes inter-site correlations. 

Without loss of variational freedom, we can impose 
that / dx,{^o\Vf 1*0) = 1 and / dx^{'^o\n^Vi |*o) = n. 
These two constraints are fulfilled by Pq = d -I- | , Pi = 
1 — 2d, and P2 = d — 6/2, having introduced the param- 
eter d = {Pq + P2)/2 and the deviation from half-filling 
S = 1 — n. The variational approach amounts to compute 
the expectation value of the Hamiltonian (1) over the state 
(2). In infinite-coordination lattices such as our Bethe lat- 
tice the computation becomes straightforward thanks to 
the parameterization in eq. (3) [16, 17], and leads to a 
variational energy per site 
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Here ho{x) = a;o/2 (—9^ + x^), {. . indicates the aver- 
age over (j)„{x), \e\ is the Fermi-sea average of the hop- 
ping energy per site. The hopping renormalization, S, is 
proportional to the overlap between the phonon wavefunc- 
tions corresponding to two different electronic configura- 
tions according to 



-S* = ^ VPTPhT / dx (l)l{x)(t)„+i{x). 



(5) 



Minimization with respect to d leads to 

U + [{hoix))o + {hoix))2 - 2(/io(x))i 

+V2au;o{{x)o^{xh)-^^^0, (6) 

while that with respect to the phonon wavcfunctions yields 
the following non-linear second-order differential equa- 
tions 

' ^ hQ{x)(f>o + V2aujoX(j)o - -^^S^/^(j)i, (7) 
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where the e^'s are Lagrange multipliers introduced to en- 
force the normalization conditions of each (f)^(x). These 
equations describe a system of forced harmonic oscillators 
coupled together by a term proportional to |e|. The lat- 
ter can be neglected in the antiadiabatic regime luq ^ |e|, 
so that the wavefunctions ^j^'s arc just displaced oscil- 
lators whose displacement is associated to the different 
local electronic configurations. In the opposite adiabatic 
regime, uo ^ |e|, the coupling term is instead dominant. 
This introduces a tight entanglement between the different 
phonon-wavefunctions and strong anharmonicity. 

As we mentioned we consider here the U s- oo limit, 
where double occupation is not allowed, hence P2 = 
and Pq = S, Pi = 1 — S, leaving us with the two coupled 
equations (7) and (8) for (j)o and (f>i. We solve them iter- 
atively by expanding the phonon wavefunctions in eigen- 
states of ho, \n), with energy nujQ as 0o = X^^nl'T-) and 
01 = J2dn\n). The equations for the coefficients c„ and 
d„ are 




Fig. 2: (Color online) m*(A)/m*(0) as a function of 7 at A 
1.5 and S = 0.3 in DMFT (squares) and GPW (circles). 
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c„' = 0, 

dn' = 0. 



In practice n ~ 100 is sufficient to get accurate values for 
the ground state corresponding to the lowest eo and ei. 

A key quantity that characterizes the electronic proper- 
ties is the effective mass m*, obtained through m* /m ~ 
1 — dYi{uj)/duj, being the single-particle self-energy 

which, for infinite coordination, is momentum indepen- 
dent. Within the Gutzwiller approach, the effective mass 
is commonly identified as the hopping renormalization fac- 
tor, which, through cq.(5), is given by 



m*(A)=m*(O)|(0i|0o) 



(10) 



where the overlap (0i|0o) includes the effect of phonons, 
while m*(0) only the correlation effects. 




Fig. 1: (Color online) Effective mass as a function of A and for 
different values of 5 at 7 = 0.2 . 



In fig. 1 we draw m*(A)/m*(0), eq. (10), for differ- 
ent dopings as function of A for 7 = 0.2, in the adiabatic 
range. The logarithmic scale on the vertical axis is used to 
show how the method reliably works for arbitrary value of 
the coupling. The linear scale in the inset emphasizes the 
onset of polaronic behavior, identified by a rapid growth 
of m*. As we anticipated, in our calculation the forma- 
tion of polarons occurs via a crossover, in agreement with 
DMFT calculations. Approaching half-filling the crossover 
becomes sharper and it moves to larger A's due to the in- 
crease of correlation effects. 

To assess quantitatively the accuracy of the GPW, in 
fig. 2 the evolution of m* as a function of 7 for A = 1.5 and 
5 = 0.3 is compared with DMFT. Like in DMFT, the vari- 
ational m* first increases with 7 in the adiabatic regime, 
reaches a maximum and finally decreases in the antiadia- 
batic regime. The GPW results closely follow DMFT for 
large 7's, and deviate more appreciably only in the adia- 
batic region, although the difference with the exact results 
is never larger than 20%. Yet, the qualitative behavior is 
correct for any 7; in particular, as previously mentioned, 
the polaron-formation remains a crossover also deep inside 
the adiabatic regime. 

The effect of the e-ph coupling on the mass renormal- 
ization is hidden in the variational calculation into the 
overlap {(f>i\(f>Q). In fig. 3 we plot |((/)i|(/)o)P as a function 
of A for two values of 7 = 0.2, 0.6 within the adiabatic 
region, where the quantitative agreement with DMFT is 
poorer. Here we compare with the same quantity from 
DMFT, and with m*(0)/m*(A) calculated in DMFT from 
the derivative of the self-energy. Obviously within DMFT 
eq. (10) has no reason to hold. In both cases the varia- 
tional wavefunction is able to reproduce extremely well the 
exact value of the overlap. What is clearly different in the 
two cases is the ability of GPW to reproduce the results 




Fig. 3: (Color online) |((;/>o|<^i)P as a function of A at 5 = 0.3 
for 7 = 0.2 and 7 = 0.6 in GPW (solid black line) and DMFT 
(dashed red line), compared with m*(0)/m*(A) from DMFT 
(dot-dashed blue line) 

for the effective mass. While for 7 = 0.6 the agreement is 
very good also as far as m*(0)/TO*(A) is concerned, which 
also means that eq. (10) is almost valid within DMFT, 
for the smaller value of 7 m*(0)/m*(A) and |((/>i|(/)o)P are 
quite different, and the accurate value of the overlap is 
not translated into an equally accurate estimate of the 
effective mass. This result is consistent with the spirit of 
our variational approach which is not in principle accurate 
for a quantity related to excitations, such as m*, while it 
correctly reproduces static quantities such as the overlap. 
Only when the phonon contribution to the mass renor- 
malization m*(A)/TO*(0) can be expressed in terms of the 
overlap, the method becomes very accurate also for this 
quantity. 

Finally, let us discuss in some details the evolution of the 
phonon wavefunctions, which are easily accessible within 
the variational approach. In fig. 4 the behavior of 4>oix) 
and 4>i (x) is shown as a function of doping (5 at 7 = 0.2 and 
for three different values of A, representative of the weak 
to intermediate to strong-coupling crossover. At A = 5.6, 
deep into the polaronic regime, the phonon wavefunctions 
are basically harmonic oscillators with a displacement very 
close to the atomic limit values (marked by the vertical 
arrows in the figure) , and with a negligible doping depen- 
dence. Similarly, when the e-ph coupling is relatively small 
(A ~ 0.8), the wavefunctions are still similar to displaced 
harmonic oscillators, although the doping dependence is 
more pronounced. The most interesting behavior occurs 
for A = 2.4, which lies in the crossover region of fig. 3. 
Here the shape of the phonon wavefunction strongly de- 
viates from the harmonic oscillator and the doping de- 
pendence is strong. In particular, 4'o{x) develops a shoul- 
der besides the main peak. The position of the shoulder 
roughly corresponds to the maximum of 0i {x) . It is ex- 
actly this deformation that makes the continuous evolu- 
tion as a function of A possible. We notice that a simi- 



lar behavior for the evolution of the phonon wavefunction 
is obtained by means of an exact numerical solution of 
the polaron problem on a two-site lattice in the Holstein 
model [18]. 
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Fig. 4: (Color online) Evolution of wavefunctions as a function 
of doping for different values of A at 7 = 0.2. The arrows 
indicate the displacements in the atomic limit. 

In this paper we have introduced a generalization of the 
Gutzwiller variational approach to account for the e-ph 
coupling besides strong electronic correlations. Wc have 
bcnchmarkcd the variational wavefunction in the infinite- 
U limit of the Hubbard-Holstcin model on the Bcthe lat- 
tice, where the variational problem can be solved analyt- 
ically and directly compared with exact DMFT calcula- 
tions. We have found that the variational results are in 
perfect qualitative, and to a large extent also quantitative, 
agreement with the exact ones. We emphasize that the 
variational wavefunction (2) treats on equal footing both 
electronic and phononic coordinates, without assuming ei- 
ther a Born-Oppcnheimcr scheme or an anti-adiabatic one. 
From this viewpoint, it may provide interesting insights 
about the way retardation effects due to the lattice dy- 
namics show up in the ground state wavefunction, espe- 
cially for strongly correlated systems where the Gutzwiller 
approach is expected to work better. 
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